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We develop a framework in which the activity of nonlinear pulse-coupled oscillators is posed 
within the renewal theory. In this approach, the evolution of inter-event density allows for a self- 
consistent calculation that determines the asynchronous state and its stability. This framework, 
can readily be extended to the analysis of systems with more state variables. To exhibit this, we 
study a nonlinear pulse-coupled system, where couplings are dynamic and activity dependent. We 
investigate stability of this system and we show it undergoes a super-critical Hopf bifurcation to 
collective synchronization. 


The collective dynamics of pulse-coupled networks of 
nonlinear oscillators has been studied extensively [T]. In 
the presence of noise, the approach has been to analyze 
the Fokker-Planck equation of the state variables of the 
oscillators. This has been highly successful for systems 
with units described by a single variable. However, it has 
proven to be difficult to extend this approach to systems 
with several state variables. This is because in such sys¬ 
tems the Fokker-Planck equation may have a highly non¬ 
trivial boundary conditions [2]. An other long standing 
theoretical framework to study irregularly pulsing units 
is the theory of the stochastic point processes. In this 
theory, the event times are described by probability den¬ 
sity functions which are history dependent. Solutions of 
the first passage time problem have long been used to 
connect this phenomenological description to the under¬ 
lying dynamics of the state variables [3]. 

In this letter, we marry these two approaches, exploit¬ 
ing the fact in that pulse-coupled systems the recurrent 
inputs into the units is fully determined by the timing 
of events. The only element that needs to be added to 
the first passage time description, is the self-consistency 
of the interactions and the network output. We first 
demonstrate our method on a simple system and show 
that the description of asynchronous state and its stabil¬ 
ity is consistent with previously derived results [3} . We 

then add a dynamic component to the interactions. This 
addition is not readily incorporated within the Fokker- 
Plank approach, however, it is easily incorporated in our 
formalism. Such dynamic recurrent couplings can be ob¬ 
served in many physical systems. For instance, tempo¬ 
ral dynamics of intracellular signaling activities is tightly 
regulated by positive or negative feedback [B], similarly 
biochemical processes concerning transmitter production 
and release in synapses in the network of neurons are 
known to be modulated by the activity of interacting 
cells 0 |B] ■ 

We consider a network of N identical oscillators with 
all-to-all feedback coupling, which receives a noisy exter¬ 


nal input. We assume that the oscillators are modeled 
as integrate and fire neurons, where their the membrane 
voltage is the state variable. Between events the (nor¬ 
malized) voltage Xi, of oscillator i satisfies 

Tm^Xi = Li{t) - Xi, ( 1 ) 

where r m is the membrane time constant and q is the 
input current into oscillator i. When the voltage reaches 
the threshold, a: t hr = 1, the oscillator emits a pulse 
and the voltage is immediately reset to the resting po¬ 
tential, x T = 0. The input, Li, can be written as 
H = G,ext + ti,fb where q jext and are the external 
and feedback input respectively. The external current is 
given by i iiext (t) = Text + er Ri{t), where the rus are inde¬ 
pendent Gaussian white noise variables, (?^(t)) = 0 and 
{ViitfVjit')) = Sij(t-t'). When a oscillator emits a pulse 
at time tk, this causes a, so-called synaptic, current input 
s(t — tk) in all oscillators. This input is given by 

s(t) = f —( e ~ t/T - - e-*M 0(f), ( 2 ) 

-/V T s i T s 2 V / 

where 0 is the Heaviside function. Here r s 2 and r s i are 
the synaptic rise and decay times. We study the network 
in the thermodynamic limit (N —► 00 ). In this limit 
the total recurrent input into all oscillators is identical, 
ti.fb = Tib, and is given by 

( 1 + Tsl Jt ) ( 1+Ts2 ^t) =9sr{t), (3) 

where r(t) = N _1 (f — f?:,fc) is the population 

firing rate. Here, ti t k is the time of the /cth event of 
oscillator i. 

To calculate inter-event density, we use that when os¬ 
cillator i emits a pulse at time t' , we have that Xi{t') = 
0 and Xi satisfies the stochastic differential equation 
T m^j.Xi = /i(t) — Xi + crriiit) until Xi reaches 1. Aver¬ 
aging over the realizations of the noise the probability 
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density p{x,t\t') for Xi(t) = x and no event has occurred 
between t' and t satisfies the Fokker-Planck equation 

= -~\}i{t)-x)p{x,t\lJ) + ^^p{x,t\lJ), 

(4) 

with initial condition p(x, t'\t') = <5(a;) and boundary con¬ 
dition p(l,t\t') = 0. The difficulty in solving this is in 
satisfying the boundary condition. For an unrestricted 
process, which satisfies Eqn. @ it is straightforward to 
show that with initial condition Xi(t') = x' the probabil¬ 
ity density p(x, t\x', t') for Xi(t) = x satisfies, for t>t' 


p(x,t\x',t') 


1 f-(x - x ave (t\x',t')) 2 \ 

(t-P) eXP V 2 J’ 

(5) 


where the noise-averaged of x is denoted as a: ave and it 
satisfies 


x ave (t\x', t') = + — [ dt" p(t")e-( t - i "V T ™ 

An Jf 

( 6 ) 


and the variance Xr is given by 
ji 


^(t) = 


<7 

2t„ 


1 — e 


— 2t/r Ir 


(7) 


The problem with initial condition p(x,t'\t') = d(x) and 
an absorbing boundary at x = 1 can be viewed as an 
unrestricted process, where a particle is inserted at x = 0 
at time t' and extracted at x = 1 at time t with some 
probability density Pr(£|P) so that 


p(x,t\t') = p(x,t\0,t’) 


dt"¥r(t"\t’)p(x,t\l 1 t"). (8) 


The inter-event probability density Pr(f|f') is determined 
by the boundary condition, p(l,t\t') = 0 and thus satis¬ 
fies the Volterra integral equation 

p(l, ^|0, t') = J dt" p(l,t\l,t")¥x{t"\t'). (9) 


The right hand side of this equation is now a convolution 
and this can be solved using the Laplace transform. The 
Laplace transform Pr eq L of Pr eq satisfies 


Ifieq.L (s) 


Peq,L (1; S|0) 
Peq,L(lj s|l) 


( 11 ) 


where /5 eq ,L is the Laplace transform of p eq . In the 
supplementary material [5], we show that the Laplace 
transform of Ornstein-Uhlenbeck density /0 eq ,L, can be 
formally calculated (for an alternative derivation re¬ 
fer to [10] ). Now, we can close the system, since 
the average inter-event interval, ( t ), can be found us¬ 
ing ( t ) = — lim s _ >0 'iP r eq,L(s)/ds and this allows to ex¬ 
press r eq = 1 / (t) in terms of p eq and <r, together with 
Peq = Pext + 5s r eq that determines r eq . 

To determine the stability of the asynchronous state 
we consider the evolution of small perturbations around 
the equilibrium firing rate, r(t ) = r eq + eri(t). With such 
a perturbation the noise averaged input into oscillators 
p{t) satisfies p{t) = p eq + epi(t), where p± is given by 


1 r°° 

Pi (f) = g s - / dt 1 ri ( t-r d -t ')(e~ 4 /Tal —e _t /Ts2 ). 

Tsl 7s2 Jo 

(!2) 

The probability density p for the unconstrained diffusion 
still satisfies Eqn. ([5]), but now with x(t\x',t') = x eq {t — 
t'\x') + exi(t\t’), where Xi satisfies 


Xi(t\t') = J dt"pi(t")t 




(13) 


Thus, we can expand p as p(x,t\x' ,t') = p eq (x,t — t'\x') + 
epi(x,t\x r ,f) + 0(e 2 ), where 

Pi{x,t\x' ,t') = -xx{t\t')-^~p{x,t- t'\x'). (14) 

ox 

Next, we write for the inter-events probability density 
Pr(f|f') = Pr eq (f — t’) + ePri(t|t') + 0(e 2 ) and insert 
this with the expansion for p in Eqn. (|9|, to obtain for 
Pri(f|t') the Volterra integral equation 


In the asynchronous state the emission rate is constant, 
r(t) = r eq , and we have p{t) = p eq = p ext + g s r eq . Con¬ 
sequently, p is given by p(x, t\x', t') = p eq (x,t — t'\x'), 
where p eq {x, t\x') satisfies Eqn. ([5| with x me (t\x',t') = 
Zave.eq^-l'l^') = p eq + (x' - p eq ) exp(- [t - t'] / T m ) . Note 
that for t —> oo, p eq (x, t\x') —t Poo(x) = exp(—(x — 
/i e q) 2 /2E^ 0 )/v / 27rE 00 , where T, x = o/y/2r m , Addition¬ 
ally, the inter-event probability density can be written as 
Pr(f|f') = Pr eq (f — t') as this density is time invariant in 
the stationary asynchronous regime. This density must 
satisfy the following Volterra integral equation 

Peq(1, 110) = [ dt'p eq (l,t-t'\l)Pv eq (t'). (10) 

Jo 


p 1 (l,t|0,t / ) - J dt" pi(l,t\0,t")Px eq (t" -t') 

= J^dt" p eq (l,t-t"\0)Px 1 (t"\t'). (15) 

Finally, we close the system using r(t), that is 

r eq +en(t) = f dt' [Pr eq (f|T)+ePr 1 (t|t , )][r e q+er 1 (t , )]+0(e 2 ). 

J — OO 

(16) 

For sufficiently small e, we can ignore terms of order e 2 . 

In this linearized system, we can make the usual Ansatz 
that r\{t) = r\e xt . With this Ansatz we will, as we 
will see below, have p\{t) = r\e xt p \, pi(x,t\x' ,t') = 
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r x e xt p x (x,t — t!\x') and Pri(t|f') = r x e xt Pr x (t — t'). In¬ 
serting this in Eqn. ([ 16 } we obtain the eigenvalue equa¬ 
tion 


poo poo 

/ dt Pr eq (f)e _At + r eq / dtPv x (t) = l. 

Jo Jo 

With r\{t) = r x e~ xt we have 


Mi = r\e 


— X t 


g s e 


— Ard 


(1 + At s i)(1 + At S 2 ) 


(17) 


(18) 


we will demonstrate that the stability analysis of these 
systems can be treated effortlessly. Here, without loss 
of generality, we focus on the couplings with depressive 
interactions (e.g. negative feedback). The network is as 
before, except that the recurrent input due to a event of 
oscillator i depends on the strength factor p it denoted as 
release factor. For instance, in network of neurons with 
depressive couplings the biophysical meaning of pi is the 
amount of vesicles that are available in synapses. The 
recurrent input is given by 


and 


xi(t\x',t') = r x e xt A x (t,t'), 


(19) 


where A At f) - ThlJS WP 

wneie A x {t,t ) — (l+Ar s i)(l+AT a2 )(l+Ar m ) • lnus ; we 

can write pi as pi(x, t\x', t') = r x e~ xt p x (x,t — t'\x') with 
P\(x,t-t'\x') = -A X (t,t')-^Peq{x,t - t'\x') (20) 


We insert this into Eqn. (15), we find hat Piq is given 
by Pri(t|t') = r x e~ xt Pr x (t — t'), where Pr A satisfies 


(i + T-si JJj (^ + T&JJj p,fb(t) = g a r r {t) (25) 

where r r (t) is the rate of release rather then the event 
emission rate. The release rate is given by r r [t) = 
uN -1 JU k pi(t~)S(t — tik), where u is the release frac¬ 
tion (see m Between events the vesicles are replenished 
with a time constant Tp, dpi/dt = (1 — p,;)rp and at the 
time of the event an amount upi of vesicles is released 
and pi is reset to (1 — u)pi. It is straightforward to show 
that if the event density is Pr(t\t'), the release rate ry 
satisfies 


dt'p eq (l,t - t'\l)e A(t * ) Pr A (t') 


= ,5a(M|0)- / dt'p x (l,t~ t'\l)Pr eq (t'). (21) 


ry(t) = u J dt'Pr(t\t')(l — e )r(f') 

+ (1 — u) f dt' Pv{t\t')e~^^^\ r (t'). (26) 


Multiplying both sides by e st and integrating over t we 
find for the Laplace transform Pr A L of Pr A 

p, ^ __ PA,l(M| 0) -PA,L(l,s|l)Preq,L(s) 

. \ 11 \ • 
Peq.LtU^T A| 1) 

Here /5a,l, the Laplace transform of p x that satisfies 


d 

P X ,-l(x,s\x ) = - — B x (p 0 ^(x,s\x‘ )-p 0:L (x, s+A+T r 


The rate in the asynchronous solution can be determined 
as follows: Given p eq and cr we calculate the Laplace 
transform Pi'l of the inter-event distribution as before. 
This determines the equilibrium rate. Using this and 


rV)) 


where B\ = 


p— At d 


(23) 

We can rewrite the 


( I+Atri ) (1 + At b2 )( 1 +Ar m ) 

eigenvalue equation (Eqn. 17) as Pro,L(A)+r eq PrA,L(0) = 


1, and plugging it in the expressions for the Laplace 
transforms obtained above, we find that the straightfor¬ 
ward eigenvalues, As, of the system that satisfy 


Eqn. (26), then we obtain for the steady state release 
rate r r , eq = ur eq /[ 1 — (1 — u)Pr L (l /t d )\. The self- 
consistency requirement is that this should agree with 
Meq = Pext+g s i"r.eq- Now, it is also straightforward to ex¬ 
tend the stability analysis to this model. We starts with 
the Ansatz that the release rate satisfies r r (t ) = r r ^ eq + 
e?7,Ae At - Following the steps of the model with static cou¬ 
plings, we obtain Pr(f|f') = Pr eq (f — t ') + ee At Pi'A(t|t') 
and r(t) = r eq + ee xt r x , with Pr A and r x proportional to 
r rj a- Combining this with Eqn. [26} and requiring self- 
consistency leads to the following eigenvalue equation 


g s r eq / dtK 1 D x (t) = K 2 / dte xt D 1 (t), (24) 

Jo Jo 

where = 1 — e“^ A+Tm « 2 = e Ard (l + Ar s i)(l + 
Ar s2 )(l + Ar m ), D x (t) = ^ (p eq (x, t|0) — p eq (x, t\l)\ x _ l 
andZ?i(t) = p eq (l, f|0)—/5 eq (l, t|l). This expression is ex¬ 
act and it directly corresponds to the eigenvalue equation 
that has been formally derived using the perturbation of 
Fokker-Plank operator j4]. 

Unlike the latter, our approach is easily extended to 
networks in which the recurrent connections are medi¬ 
ated through coupling with dynamic strength. In this part 


g s r eq ^u J dtKiK 3 D x (t) + (1 - u) J dt n^D^tJj = 

/ poo poo 

J dt K 3 e~ xt Di(t) + (1 — u) J dt K 4 e~ xt Di(t) 


where, n 3 = \ — e~ t l TD and K 4 = e.~ t l TD . Now, we can as¬ 
sess the stability of the network with dynamic couplings, 
by finding As that satisfy the above equation. We numer¬ 
ically determined the eigenvalues for different Tp and u, 
adjusting g s to jeep the rate constant in the asynchronous 
state. The asynchronous state destabilizes through a pair 
of purely imaginary As, corresponding to the emergence 
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FIG. 1. Phase diagram of a system with depressive couplings. 
The asynchronous stationary state is stable underneath the 
curve. The solid line correspond to the parameter regime 
where the purely imagery eigenvalues give rise to Hopf bifur¬ 
cation of asynchronous irregular state. The marked symbols 
are the parameters in the phase space that we adopt to numer¬ 
ically evolve the system [§] in Fig.2. Parameters: p ex t= 0.95, 
er=0.0228, g s =0. 00245, r m =0.020 and t S 2=0.001, a: t hr=l and 
14=0. 
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FIG. 2. The numerical simulation of the full system with 
depressive couplings using population density treatment [9]. 
Each subplot illustrates the population firing rate of the sys¬ 
tem with the parameters are marked with corresponding sym¬ 
bol in Fig.l. 


of s limit cycle oscillations due to an Andronov-Hopf bi¬ 
furcation. Figure 1 shows the resulting phase diagram. 

The approach introduced here can be utilized to evolve 
a system for N —> oo, as the time dependent inter-event 
density can be self-consistently determined. By exploit¬ 
ing this, we numerically evolve a network with depressive 
coupling [S] to study the behavior near the bifurcation 
point. The activity continuously changes after the bi¬ 
furcation point (Fig.2) and the amplitude of collective 
synchrony grows indicting a super-critical Hopf bifurca¬ 
tion [9|. 

In the present letter, we derived self-consistent descrip¬ 
tion of the inter-event distribution of non-linear pulse 


coupled oscillators with interactions. We additionally 
characterized the asynchronous state and its stability. 
For static interactions, this result coincides with the re¬ 
sult from the classical Fokker-Planck approach. How¬ 
ever, as we showed our approach is easily extended to 
incorporate the effect of dynamical coupling. Using this, 
we investigated how networks with interaction through 
couplings with short-term depression undergoes a Hopf 
bifurcation to a state with collective synchronization, so- 
called population spikes. The method also allows for a 
efficient way to simulate the dynamics of systems in their 
thermodynamic limits [5]- We showed in this limit the 
system may exhibit a super-critical Hopf bifurcation. Up 
to now, the collective effect of activity dependent mod¬ 
ulation of interaction in a network has only been ana¬ 
lyzed in models without non-linear contributions of the 
inter-event density mm- In such networks interesting 
phenomena such a as Shilnikov chaos has been observed 
when positive feedback (e.g. facilitation) are added flTj . 
It is straightforward to extend our approach to study 
those cases where the network of oscillators with self- 
consistent inter-event density (e.g. spiking neurons) are 
also considered. The method presented in here relies on 
our ability to calculate p , the solution of the unrestricted 
of Fokker-Planck equation. Once this is achieved impos¬ 
ing the necessary boundary conditions is easy using the 
presented method. Thus we believe that our approach 
will have a wide range of applications. 
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